Molecular dynamics modeling of the effects of cementation on the acoustical 
properties of granular sedimentary rocks. 
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The incidence of cementation processes on the acoustical properties of sands is studied via molec- 
ular dynamics simulation techniques. In our simulations, we consider samples with different degrees 
of compaction and cementing materials with distinct elastic properties. The microstructure of ce- 
mented sands is taken into account while adding cement at specific locations within the pores, such 
as grain-to-grain contacts. Results show that acoustical properties of cemented sands are strongly 
dependent on the amount of cement, its relative stiffness with respect to the hosting medium, and 
its location within pores. Simulation results are in good correspondence available experimental data 
and compare favourably with some theoretical predictions for the sound velocity within a range of 
cement saturation, porosity and confining pressure. 
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I. INTRODUCTION 

After settling, mechanical and permeating proper- 
ties of sand are modified by processes called diagenesis. 
Among other diagenetic processes, cementation is denned 
in a broad sense (Bryant, 1993) as the introduction of 
mineral phases into rock pores. Many mechanisms of ce- 
mentation are discussed in the literature: minerals can 
be transported in suspension with underground waters 
and settle in the pores of a final hosting medium. In low 
flux zones, cements can form as a solid precipitate on 
grain surfaces. Formation of clay bonds between grains 
and grain interlocking as a product of the matrix dissolu- 
tion (Gundersen et al., 2002) are among others, possible 
cementation mechanisms. 

Depending on their origin and location within the 
pores, cements modify hosting rock properties in differ- 
ent manners. It has been experimentally observed that, 
if cement is deposited near grain-to-grain contacts, it can 
significantly increase strength and stiffness of a granular 
material (Bernabe et al., 1992). In this case, grain rota- 
tion and displacement are inhibited due to cement bonds 
and samples develops a frame that resists compaction. In 
an experimental study with glass beads (Hezhu, 1992), 
the presence of epoxy cement at contacts was observed to 
prevent bead crushing during compaction. Theoretically, 
this effect was attributed to a more uniform stress dis- 
tribution within contact areas due to cementation. This 
kind of cement is the main responsible of rocks' cohesion. 
Therefore, the lack of this kind of cement is a potential 
handicap for the stability of boreholes. 

A recent experimental study has addressed how ce- 
ments can modify hosting rock microstructure and re- 



lated hydraulic properties, such as permeability (Prasad, 
2003). According to the authors, microstructure can 
be characterised through a geometric index FZI = 
l/\fP\TS v , where P s is a pore shape related value, r 
is the tortuosity and S v the specific surface area per unit 
of grain volume. Such geometric index, is then correlated 
with permeability and wave velocity. As an example, if 
cement is deposited uniformly on grains surfaces, specific 
surface area sharply increases. On the other hand, pore 
filling materials can also grow from grain surfaces towards 
pores, forming intertwined structures that increase tor- 
tuosity (Manmath and Lake, 1995). Both mechanisms 
lead to different microstructures at the pore scale that 
affect rock hydraulic and mechanical properties. 

Despite the importance of the prediction of cement ef- 
fects, this is often a difficult task since the exact location 
of cement in pores is usually unknown. Although some 
insights can be obtained from microphotography or other 
high resolution methods, the usage of these techniques 
is limited when considering rocks at larger (formation) 
scales. As an alternative, acoustical methods are of great 
value in the characterisation of the underground at larger 
extents, being successfully used in the identification of 
underground structures and fractures. Bowever, for the 
precise interpretation of acoustical logs it is necessary to 
advance in the understanding of the acoustics of porous 
media and identify the connections between macroscopic 
acoustic observables, and microscopical petrophysical pa- 
rameters. In this sense, due to the apparent linkage of 
rock stiffness and cementation, it seems plausible to use 
acoustic methods to estimate the type and amount of 
cement present in rocks. For this purpose, one must es- 
tablish the relation between the acoustical observables, 
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such as velocity and attenuation, and the cement satura- 
tion, its elastic properties and localisation within pores 
among other parameters. 

Experimentally, the acoustical effects of coating ce- 
ment were studied when considering artificially fabri- 
cated samples cemented with sodium silicate (Nakagawa 
and Myer, 2001). Samples showed increasing wave ve- 
locities and attenuation with cement saturation. The 
acoustical effects of cement located near contacts, were 
studied experimentally by considering epoxy cement and 
glass beads (Hezhu, 1992). In this work, wave velocities 
sharply increased with cement saturation. Same results 
were reported in (Dvorkin et al., 1994a; Dvorkin et al., 
1994b) and references therein. 

From a theoretical point of view, Dvorkin et al. (1994a) 
proposed a model for cemented sands based on the effec- 
tive medium approach. This model considers the elas- 
tic properties of a granular material whose porosity de- 
creases by the addition of cementing material near grain 
contacts. The model is based on the analytical solution 
of the related elasticity problem. Theoretically, it is ex- 
pected that small amounts of cement added precisely at 
grain contacts, to sharply increase rock stiffness. Differ- 
ent cementation schemes, where cement is added prefer- 
entially far from contacts, or as a coating on grain sur- 
faces were considered in (Dvorkin and Nur, 1996) and the 
results summarised by Mavko et al.(1988). For reference, 
the main formulae are also presented here in appendix^ 
This theoretical formulation was proposed as an expla- 
nation for some high-velocity high-porosity samples and 
different velocity-porosity trends observed in North Sea 
sandstones (Dvorkin and Nur, 1996). 

Numerical methods, such as molecular dynamics or fi- 
nite element methods, represent a new approach to the 
study of complex multibody systems such as sand. In this 
approach, few numerical studies consider the effects of ce- 
mentation on the mechanical and hydraulic properties of 
sands. In early papers, porous media was modeled us- 
ing disks bonded by elastic springs representing cement. 
The mechanical properties of this model were observed 
to depend on the bonding distribution and properties. A 
three dimensional approach studied the brittle failure of 
cemented rocks under external load, modeling cement as 
breakable grain bonds (Guodong et al., 2003). 

Previous numerical studies considered cementation ef- 
fects on permeability by emulating the microstructure of 
cemented rocks and then solving flux equations in this 
geometry. Some approaches consider sand as a pack of 
spherical particles whose radii is extracted from a given 
distribution. The cementation process is then simulated 
by increasing grains radii by a given amount, modeling 
the quartz cement overgrowth (Bryant et al., 1993; Bakke 
and Oren, 1997). More recent studies account for prefer- 
ential growth of cement to pores or throats by a more 
elaborate algorithm that simulates cement overgrowth 
from grain boundaries to the surface of their Voronoi 
polyhedron (Schwartz et al., 1987). The resulting spa- 
tial correlations and other geometric characterisation was 



studied numerically by Biswal et al.(1999). 

However, to the best of our knowledge, acoustics of 
sands is an open problem where numerical methods can 
be very useful to understand the underlying physics be- 
hind the phenomena, to test theoretical models and to 
improve their predictive capabilities. In our approach, 
we move a step further in the modeling of acoustical re- 
sponse of cemented sands. We propose a simple particle 
based method (molecular dynamics) that considers the 
amount of cement and its location within pores. 

Particle based methods or discrete methods differ from 
other well known techniques such as finite elements in 
that no continuous wave equations are solved directly. In- 
stead, molecular dynamics mimics the underlying physics 
of wave propagation at the micro or meso scales using 
micromechanical interaction rules between discrete parti- 
cles. This method involves solving Newton's equations of 
motion for an N particulate system and global behaviour 
is obtained by the cooperative effect of interacting parti- 
cles (Rappaport, 1995). 

In our simulations we study the acoustical effects of 
soft and hard cements on a model sample. We consider 
the case when cementing material is added preferentially 
to grain-to-grain contacts, on grain surfaces or as a solid 
material added preferentially to pore bodies. Our model 
sample mimics a well sorted naturally occurring sand, 
and cementation is simulated capturing the underlying 
physics governing the acoustic pulse propagation on ce- 
mented sands. 

This paper is organised in three main parts: section 
ITT1 summarises the basic concepts related with cementa- 
tion process. The following section ITTT1 covers a detailed 
description of the methods used in the simulations. Fi- 
nally, simulation results are presented and discused in 
section IIVI Important information regarding with the 
calibration of simulations and some theoretical formulae 
are presented in appendix. 



II. CEMENTATION PROCESS 

Clastic rocks are often cemented by addition of clay 
particles. These cements can be classified in two broad 
categories: Detritical and Diagenetical (Manmath and 
Lake, 1995). The former are transported in suspension 
with subsurface waters to a final location in a hosting 
medium where they deposit. This kind of cementing 
particles are mainly attracted to the sharpest zones of 
pore space, such as grain-grain contacts, where the sur- 
face forces are largest (Adamson and Gast, 1982). These 
cements reinforce grain contacts, increases rock's stiff- 
ness and sound propagation velocity. Theoretically is 
expected that even small amounts of soft cement added 
precisely at contacts, to sharply increase wave velocity 
propagation (Dvorkin et al., 1994b). 

On the other hand, diagenetic cements form as precip- 
itates of solid material initially in solution in meteoric 
waters (or low mobility waters) . Due to its chemical ori- 
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gin, many varieties of these clays are possible, depending 
on the concentration of different species in pore fluids 
and chemical equilibrium at underground pressure and 
temperature. In this category, the so called pore lining 
clays deposit preferentially as coatings on grain surfaces, 
except at the grain to grain contacts (Guodong et al., 
2003). The most frequently occurring pore lining cements 
are smectite, illite and chlorite. This kind of cements, 
if deposited homogeneously on grain surfaces, leads to 
a sharp increase in specific surface area of pores while 
porosity reduces. Underground sections with this kind of 
cementation can exhibit a decrease in permeability under 
the normal compaction trend (Prasad, 2003). 

Another possibility for clays is to grow as long crystals 
from grain surfaces to pores. This kind of clays, such 
as radial illite and smectite, are sometimes called pore 
bridging. The effect of these cements is mainly to re- 
duce porosity and permeability while increasing tortuos- 
ity. As cement is mainly deposited away from contacts, it 
does not significantly contribute to sample stiffness until 
porosity is sufficiently reduced. 

Other cementing materials, such as quartz or calcite, 
are also present in rocks as cements. These cements can 
locate preferentially at pore bodies (pore filling cements) 
forming tight structures. Quartz overgrowth is the for- 
mation of quartz crystals on the sand grains surface. Al- 
though crystals exhibit euhedral faces, the overgrowth 
approximates also to concentric rims or coatings. 

III. THE MODEL 
A. Settling Process 

The initial grain pack is built by following a ballis- 
tic algorithm to simulate the settling process of clastic 
sands (Bakke, 1997). The procedure models sand grains 
as spherical particles that fall on an initially flat surface 
from a fixed height in the z direction and random co- 
ordinates in the xy plane. The radius of each particle 
is extracted from a given grain size distribution. In our 
simulations, the radius R of each particle is selected with 
equal probability in the range 0.018 — 0.02cm, modeling 
a well sorted sand. 

During settling, Newton's equations of motion are 
solved for grains falling under the influence of a grav- 
ity field ~gz, a viscous drag due to the aqueous medium 
and the interaction with other previously settled grains. 
At this stage, when two grains come into contact, they 
interact with a repulsive non-linear viscoelastic force F c 
given by: 

F c = {(/<„£ - 7nC}n + F s . (1) 

The first term in the Eq. represents the force F n 
normal to the contact area. For the labeled grains 1 and 
2 with radii R\, R2 and positions ri, r2, £ = max(0, R± + 
R2— |i"2 — ri|) and n.12 = ri — rg/lri— r2 1 is the unitary vec- 
tor joining the grains centers (see Figure^. The normal 



stiffness of the contact K n = i/3-J£RfEf is computed ac- 
cording with Hertz theory, where Rf = (R1R2) / (R1+R2) 
is an effective radius and a c — W£Rf is the radius of con- 
tact area. The effective Young's modulus Ef is calculated 
as 1/E f = (1 - vl)/E x + (1 - 4)/E 2 , being E u E 2 , v x 
and V2 the Young modulus and Poisson's coefficient of 
the grains material (Love, 1944). The second term in 
braces in Eq. ^ represents the the viscous forces for nor- 
mal deformations, where j n = JniS.Rf) 1 ^ 2 and 7„ is a 
damping constant. 




FIG. 1: Interaction forces are decomposed in normal and 
tangential directions to the contact area. 



F s = -min (\k s (\, fJ,\F n \) sign(Qs, (2) 

The shear force F s in Eq. ^ depends on contact his- 
tory and cannot be entirely determined by the position 
of grains. Such force is given by Eq. where the tan- 
gential stiffness is n s = 8a c G/(2 — is), shear modulus of 
the grains is given by G and v is the Poisson's coeffi- 
cient. The term a c denotes the radius of contact area 
previously defined and fi represents the Coulomb friction 
coefficient. The term £ denotes the component of the 
relative displacement vector in the tangential direction 
that took place since the time to, when the contact was 
established. 

at) = / v s (t')df, (3) 

J to 

where v s = v.y ■ s and s = C/lCl- The vector £ can be 
computed from the displacement vector d between the 
grains £ = d — (d • n )n (see Figure^]). 

At this stage, periodic boundary conditions are im- 
posed on gravity perpendicular directions to emulate an 
effectively larger system and the underground natural 
confinement of rocks in those directions. 

The initial pack models a high porosity cf) — 41.2% 
sand, where contacts involve marginal contact areas. The 
size of the generated sample was 8 x 8 x 32 in units of max- 
imum grain diameter. Following previous results (Garcia 
et al., 2004), this volume is several times larger than the 
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minimum homogenisation volume for the porosity of a 
sample with the given grain size distribution. 

Samples with different packing are obtained by simu- 
lating the naturally occurring compaction of sediments. 
Following Garcia et al.(2004), this process is simulated 
by displacing towards the sample core, grains contained 
in frozen slabs at both ends of the sample in the z direc- 
tion. Once such macroscopic strain in imposed, ample 
time is allowed for grains to rearrange and relax accumu- 
lated stress. 



B. Cementation model 

Once the initial sand pack is obtained, the cementation 
process is simulated in several steps representing different 
stages in the sand diagenetic history. At each step, a vol- 
ume of cement is added at specific locations within pores. 
Cement is added in three different manners that model 
cement with distinct origins and their micro-geometric 
implications. 

The procedure followed to add cement starts by dis- 
crediting the three dimensional sample into cubic boxes 
of side length c w ~ i?/10, R being the average radius of 
matrix grains. The set of cells intersects the sample voids 
and solid matrix, but just those cells whose centre coin- 
cides with pore space are potential locations for cement. 
Cells are subdivided into three groups: 

1. Contact cells are those that intersect two or more 
grains. These, are located near grain contacts at 
the sharpest zones of pore space. 

2. Surface cells are those centred in pore space and 
contacting with at least with single grain. Contact 
cells are a subset of surface cells. 

3. Body cells are those centred in the pores without 
intersecting any grain. 

In practise, cement is added by marking as cement 
saturated a fraction of pore centred cells. This mark is 
achieved by locating a cement particle at the centre of the 
cell. Each cement particle represents a cubic subvolume 
in the sample where cement is allocated. 

C. Contact Cement 

To model the microgeometry associated with contact 
cement model, we propose an algorithm of several steps 
that represent a different stage in the diagenetic history 
of the sample. Each step involves porosity reduction due 
to the addition of cement to pores. First, a small fraction 
of contact cells is marked as cement saturated. During 
the following steps, an increasing number of contact cells 
are filled with cement until no contact cells are left. 

In practise, a cell are marked as cement saturated 
by adding a single cement particle to the centre of the 



cell. As this particle represents a cubic subvolume 
of cement within the sample, we assign to it a mass 
m c = p c c%, where p c is the volumetric density of the 
cementing material. 

After all contact cells are filled, further cementation 
steps involves the addition of cement to selected cells, 
neighboring those previously filled. First, we list all 
empty neighboring cells of those already filled with ce- 
ment. Then, among these empty cells we select those 
lying closer to the grain matrix, and fill them with a new 
cement particle. Thus, when all contact cells are filled, 
this procedure promotes the cement accumulation from 
contacts to surface cells. 

A sequence of steps where cement was added according 
to the described procedure is shown in Figure the 
Figure shows a two dimensional example for illustrative 
purposes. However, results reported here refer to three 
dimensional samples and smaller cells (cells width c w ~ 
1 th of average grain radius) . 




FIG. 2: Figures A,B: an increasing number of contact cells 
is filled with cement until no contact cells left. Then, ce- 
ment grows from contacts to surfaces in Figures C to D. In 
the Figure, 2D samples are shown for clarity. Nev- 
ertheless, results reported here refer to the three- 
dimensional samples. 



D. Coating cement 

To model the microstructure related with coating ce- 
ment, we mark each cell with an index proportional to its 
distance to the solid surface. At each cementation step, 
all empty cells are checked and those with the smallest 
index (closest to matrix surface) are filled with a cement 
particle. In the first steps, cementing material is added 
on grain surfaces. As is shown in Figure further ce- 
mentation steps add cement deeper into the pore bodies 
so that rims appear around grain surfaces. 

The successive cementation steps simulate several 
stages of the cementation history of the sample, where 



5 




FIG. 3: From step a to c cement is added to grain surfaces 
as a coating. This procedure simulates the homogeneous pre- 
cipitation of solid material on the surface of pores. Figure 
shows a two dimensional case for clarity, reported re- 
sults are for three dimensions though. 



a change in the chemical equilibrium conditions causes 
precipitates to deposit uniformly on the grain surface. 

E. Pore body cement 

Simulation of pore body cements begins with the ad- 
dition of cement particles to a fraction of surface cells 
selected at random. In the following steps, we check 
the empty cells contacting those already filled with ce- 
ment. Each empty cell is then filled with a cement parti- 
cle according to a probability p, proportional to the cell 
distance to grain surface. This procedure, simulates a 
preferential grow of cement from grains surface to pore 
bodies. 

As shown in Figure^] the initial cemented cells behave 
as seeds for clusters that grow from the grain surface 
to pore bodies. The described procedure guarantees the 
continuity of the cement phase and simulates cementing 
materials in the form of long crystals attached to grain 
surfaces. 




FIG. 4: Figure a: initial seed of cement. Figures b,c: fur- 
ther cementation steps when cement is added preferentially at 
pore bodies. Figure shows a two dimensional cases for 
clarity. However, results refer to three dimensional 
samples. 

Cell width defines the resolution limit of the simula- 
tions. If cells are too wide, poor resolution is obtained 
and porosity decreases sharply to zero after a few ce- 
mentation steps. Small cells allow more control of the 
final porosity of the sample at expenses of an increase 
in the computational cost. More details about the effect 
of resolution on the computational cost are found in the 
appendix [5] 



F. Interaction 

After a given porosity is reached due to cementation, 
we simulate the acoustical pulse transmission through the 
sample using Molecular Dynamics techniques. As Molec- 
ular Dynamics involves the explicit solution of Newton's 
equations of motion for a particle based system, we need 
to consider the interaction forces between the particles of 
the different constituents of the aggregate. In our model, 
particles conform the quartz matrix and cement. We 
do not consider fluid particles that may be in pore vol- 
ume. This approximation is justified if we suppose the 
fluid saturating pores to be air free to drain while sample 
strains. The effect of a more incompressible saturating 
fluid such as water or oil could be included and it is 
proposed for future work. In this way, our simulation 
techniques allow us to isolate the acoustical effect of ce- 
ment from other agents that modify the acoustics of a 
real sample. With this approximation, we need only con- 
sider three kinds of interaction forces between particles: 
grain to grain forces, interaction between cement parti- 
cles, and interaction forces between cement particles and 
grains. 

In our model, the addition of cementing material is 
assumed not to disrupt the internal stress state of the 
sample. This assumption is justified since cementation is 
a relatively slow process when compared to compaction 
(Bernabe et al., 1992). Therefore, cement is added to 
pores at mechanically stable locations after grain-to- 
grain contacts are formed. When a perturbation of this 
equilibrium state occurs, changes in the internal stress 
field of the sample are produced. At the pore scale, 
grains change their relative positions and interparticle 
forces change by an amount AF C with respect to the ini- 
tial equilibrium network. 

If after a perturbation a pair of contacting grains la- 
beled 1 and 2 at positions ri and r 2 , displace to ri + dri 
and Y2 + dr 2 respectively, their relative displacement is 
calculated as d = dr2 — dri. The component of the 
displacement along the normal direction of the contact is 
given by A£ = dn, being n = (r 2 — ri)/|r 2 — ri| the unit 
vector joining grains centre at equilibrium. The tangen- 
tial displacement is calculated as A( = d — (d- n)n. This 
vector defines the tangential direction § = A£/A£ with 
AC = |AC|. During this stage, the force change AF C due 
to relative displacement of contacting grains is computed 
according to: 

AF c = {(fc n A£- 7 „A£)}n + AF s . (4) 

Term K n in Eq. 0] is the elastic stiffness for normal de- 
formations, defined previously in equation ^ However, 
in this stage we assume a negligible change in the con- 
tact area due to the acoustic pulse. This approximation 
constitutes a linearization of Hertz force, that represents 
the normal grain— grain elastic interaction that acts on 
contact. The second term in braces in Eq. 0] describes 



TABLE I: Parameters used in Eq. [I]and[S]to compute forces 
between quartz grains (Mavko et al., 1988). 



Parameter 


Symbol Value 


Density (gr /cm 3 ) 


Pa 


2.65 


Shear modulus(GPa) 


G 


44 


Poisson coefficient 


V 


0.08 


Friction coefficient 


/< 


0.3 



the viscous forces for normal deformations. In our sim- 
ulations, we have adjusted this parameter to obtain a 
normal restitution coefficient e„ = 0.9, in agreement with 
experimental results (Kuwabara and Kono, 1987; Schafer 
et al., 1996). 

The term AF S in Eq. ^represents the change in shear 
force. To compute this term, we assume that no sliding 
between grains occurs during pulse transmission. In such 
case, shear force is considered viscoelastic and calculated 
according to 

AF S = {fc s AC- 7 sAC}§. (5) 

The first term in the Eq. represents the elastic com- 
ponent of the shear force. As no sliding is assumed, the 
tangential stiffness is computed according with Hertz- 
Mindlin theory k s — 4a c G/(2 — v). The term G denotes 
the grains' shear modulus, v is the material Poisson's ra- 
tio and a c is the radius of the initial contact area defined 
previously. The second term in the Eq. \5\ represents 
the viscous forces for shear deformations. Parameter 7 S 
was adjusted to obtain a tangential restitution coefficient 
e s = 0.9. The parameters to compute forces, were chosen 
to model quartz sand grains (see Tabic |D . 

To compute forces between cement particles, we follow 
a procedure described by O'Brien and Bean (2004). In 
this procedure, cement particles are considered as nodes 
on the cubic lattice that discretise the sample. When a 
cement particle displaces from its initial position on the 
lattice, it interacts with its 18 closest neighbours in the 
lattice as shown in Figure [SJ This interaction force is 
calculated according to: 

C 

f c = C n [ui2.fi] n + pjUi2, (6) 

I x i2| 2 

where X12 denotes the vector joining particles 1 and 2 on 
the undisturbed lattice and n = X12/IX12I. The vector 
U12 denotes the relative displacement of particles from 
their relative positions on the undisturbed lattice. The 
first term is the normal force with elastic constant C n . 
The second term in Eq. 0is the bending term and C is 
an elastic constant for the interaction (O'Brien and Bean, 
2004) . The force given in Eq. [5] can be expressed in a 




Neighbors 



FIG. 5: Each cement particle interacts with cement particles 
present in each of its closest eighteen nodes. Interaction force 
is decomposed in a component that acts along the line joining 
cells centre and a bending term. 



form similar to Eq. 0] However, we decided to remain as 
close as possible to the original reference for clarity. 

Following O'Brien and Bean (2004), energy of the elas- 
tic lattice can be compared with the elastic energy of an 
elastic continuous. The shear elastic modulus G c and 
Lame constant A c of the lattice can be expressed in terms 
of C n and C according with the relation: 

\ _ C n 2C 

C ~~ Cu, c3 > 

(7) 

r> — c„ , 2C 
c ~ ~ ^' 

where c w is the lattice parameter (in our simulations 
c w w 1 th of average grain radius) . 

V 2 — -L ( + 2C \ 

P C Pa \ C W ~^ ~cZ) ' 

(8) 

From linear elasticity, it is possible to express the wave 
velocities in the cement phase as in Eq. where V pc and 
V sc are respectively the P— wave and S— wave velocity 
of the cementing material. The term p c is the volumet- 
ric mass density of cementing material. As each cement 
particle represents a cubic subvolume of cement, such 
density is related to the mass of cement particles m c . 

C = \ Pc (W pc 2 - V sc 2 ) cj (9) 
m c = p c c w . 

For a cement of given mass density p c and wave veloc- 
ities Vpc, V sc , the mass of cement particles m c and elastic 
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TABLE II: Elastic parameters and mass density used to model 
different cementing materials. 



Parameter Symbol Hard cement Soft cement 



P- Velocity {K m/s) 




3.0 


1.5 


S-Velocity(/srm/s) 


v sc 


1 .8 


0.9 


Density (gr /cm 3 ) 


Pc 


2.65 


2.65 


Poisson Coef. 


v c 


0.22 


0.22 



constants in Eq. |S|are determined through Eq. |5J In our 

simulations, we considered the hard cement and soft ce- 
ment limits compared to quartz. Parameters used in both 
cases are shown in Tabled The hard cement case, mod- 
els a quartz-like cement such as that formed by quartz 
dissolution of matrix grains. Soft cement models a ma- 
terial added to pores such as soft clays. To simplify, we 
considered the same volumetric mass density for cement 
in both cases. 

The third kind of interaction force arises when grains 
displace relative to cement particles contacting their sur- 
face. This interaction in also computed by Eq0J but the 
elastic constants n n and k s are now computed consider- 
ing the Young's modulus and Poisson coefficient of quartz 
and cement. The diameter of contact area is taken as the 
cement cell width c w . In this manner, a grain feels a force 
due its contacting grains and due to cement particles con- 
tacting its surface. 

G. Acoustics 

In order to perform acoustic tests, the sample is subdi- 
vided along the propagation axis z into parallel slabs of 2 
largest grains diameters. A square wave displacement is 
imposed on grains contained in one end slab of the sam- 
ple while the other ending slab is kept fixed as is shown in 
Figure EJ The former boundary condition, emulates the 
setup of an acoustic measurement. The fixed boundary 
condition on the opposite end prevents the sample from 
displacing as a whole. 

The imposed strain displaces grains from their original 
equilibrium positions and disrupts the initial force net- 
work. This perturbation propagates in the form of an 
acoustic pulse with a well defined central frequency. To 
keep track of the pulse position and amplitude, we moni- 
tor changes in the forces felt by grains as function of time 
and position. The average of the force along z over the 
grains contained in each slab at time t is taken as the 
acoustical signal A(z, t) at the position of the slab centre 
z. This value is calculated according with the relation: 

1 N 



Forced S| a bs 




FIG. 6: Sample volume is discretised into parallel slabs along 
z axis. At one end, a square wave displacement is imposed 
while the opposite slab is frozen. Average forces on grains 
contained in each slab are calculated as function of time to 
track pulse position and amplitude. 



where term /j is the force change at time t on the ith grain 
or cement particle in the slab. The summation index runs 
over the N particles in the slab. This procedure simulates 
a set of acoustical detectors distributed along the sample. 

As Figure [7| shows, the position in time of the stress 
pulse can be tracked by comparing signals received in 
different detectors. 



4: 



t 




Time (jis) 

FIG. 7: Signals recorded, at the simulated acoustical detec- 
tors, correspond to the average force along z of grains and 
cement particles in different slabs of the sample. From top to 
bottom, signals are recorded at detectors located increasingly 
farther from the source. 



L LAT + TAL 

Mean pulse velocity V p is calculated according to Eq. 
1111 by comparing the arrival times T\ and Ti of the signal 
maximum in two detectors separated by a distance of 
L. The error in pulse position AL « R is estimated 
as the half width of slabs *=s 2R. Error in time arrival 
AT, is taken as the time spacing between two consecutive 
calculations of force average in the slabs. 
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IV. RESULTS 

We simulate the acoustic pulse propagation through 
a cemented sand. Cement is preferentially added near 
grain-to-grain contacts following the procedure described 
in section Till CI and its elastic properties are given in the 
column labeled as Hard Cement in Table |H] These pa- 
rameters are suited to model a quartz-like cement and 
common used cements in experiments. This initial sys- 
tem mimics a well sorted and non compacted sand with 
porosity <f) = 41.2% and point grain-to-grain contacts. 
We refer to this system as the non-compacted sample. 
Porosity reduction of the initial high porosity sample is 
achieved in several steps in which varying amounts of ce- 
ment are added within pores. Acoustic pulse propagation 
is simulated for different degrees of cement saturation S c . 

3.5 n ' 1 ' 1 ' 1 ' 1 



3 - 




0.1 0.2 0.3 0.4 

Cement Saturation S„ 



FIG. 8: Small amounts of contact cement contributes to in- 
crease stiffness of the composite. As more cement is added, 
it locates farther from contacts and stiffness increases at an 
smaller ratio. Consequently, V P (S C ) trend shows a character- 
istic downward concavity. 

As shown in Figure [SJ we observe an increase in sound 
velocity when an small amount of cement is added pref- 
erentially near grain contacts. This trend occurs be- 
cause cement increases the restoring force of contacts and 
bridges grains initially separated. As a result, stiffness of 
the sample increases at a higher rate than porosity is 
reduced. Consequently, V p trend shows a downward con- 
cavity. As has been pointed before by Dvorkin and Nur 
(1996), this is a consistent explanation to some observed 
high porosity-velocity samples. 

When more cement is added, porosity reduces and par- 
ticles tend to locate progressively farther from contacts. 
In this case, stiffness increases at an slower rate than in 
the higher porosity range, when cement particles locate 
precisely at contacts. The net result is a decreasing slope 
of the velocity-porosity curve as porosity is reduced. 

To compare our results with effective medium approx- 
imations (Dvorkin et al., 1994a; Dvorkin and Nur, 1996), 
we plot in Figure 03 the theoretical bounds for the sound 



velocity of the cemented sample. Both theoretical curves 
were obtained with the formulae given in appendix 1X1 for 
the simulation parameters. The upper trend (high ve- 
locity) corresponds to the idealised case when all cement 
is precisely added at contacts. The low velocity trend, 
corresponds to the scheme in which cement is added uni- 
formly on the grain surface. 

As shown in the Figure |SJ obtained velocity-porosity 
trends for the non-compacted sample with contact ce- 
ment, have the same qualitative behaviour than theoret- 
ical predictions. However, theory predicts a sharp ve- 
locity rise when small amount of cement is first added. 
Beyond the initial velocity rise, theoretical curves tends 
to saturate. Instead, we observe a soft continuous veloc- 
ity increase with cement saturation S c . 

For a set of parameters close to ours, previous experi- 
mental studies reported a similar qualitative behaviour 
for a system of beads cemented with solidified epoxy 
(Hczhu, 1992; Dvorkin et al., 1994a), but the experimen- 
tal results lay between both theoretical curves (Dvorkin 
et al., 1994a). We attribute this discrepancy with our 
results, to the stiffness of the initial uncemented sam- 
ple. Theoretically, a finite non vanishing velocity is ex- 
pected for the uncemented samples regardless of the ef- 
fects confining pressure. For the case shown in the Fig- 
ure [SJ the minimal theoretical velocity is higher than 
that obtained in our simulations. Experimentally, re- 
sults in the aforementioned experiments deal with com- 
pacted samples, that uncemented, propagate sound at 
velocities higher than the minimum predicted in the the- 
ory (Hezhu, 1992). On the other hand, in our settling 
algorithm grains settle due gravitational forces and no 
confining pressure is applied. As a result, we obtain very 
loose samples with marginal contact areas where P-Wave 
velocity can be quite smaller. 

With the aim to better reproduce common experimen- 
tal setups used for rocks, we simulated the compaction 
of the uncemented sample to a porosity <fi = 38.6%. At 
this porosity, confining pressure was recorded as P c ~ 
SQMPa and sound velocity was observed to increase from 
V p = 0.2Km/s to V p = l.7Km/s, similar to the mini- 
mal velocity reported in the aforementioned experiments 
(Dvorkin et al., 1994a). Then, the cementation process 
was simulated as in the previous case and acoustic tests 
were performed for different cement saturations. 

As expected, Figure shows that for zero cement sat- 
uration, velocity is higher in the compacted sample than 
in the uncompacted sample. Additionally, the Figure 
shows that in the case of the compacted sample, velocity 
increases at a slower rate when cement saturation S c in- 
creases. The net result is that the velocity trend for the 
compacted sample lies between theoretical curves. In this 
case, for the range of cement saturation S c shown in the 
Figure^ obtained results are in good quantitative agree- 
ment with available experimental data reported in pre- 
vious works (Hezhu, 1994; Dvorkin et al., 1994; Dvorkin 
and Nur 1996). It should be noted that for the case 
of the compacted sample, theoretical curves are slightly 
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different from those of the uncompacted sample, since co- 
ordination and porosity are different. Nevertheless, these 
variations are relatively small. 

These results suggest that confining pressure of the ini- 
tial uncemented sample is a crucial parameter which de- 
termines the velocity trends obtained after adding small 
amounts of cement. This is an intuitive result since for 
loose sediments close to suspension, stiffness can be neg- 
ligible. Then, any stiffness increase due to cementation 
can be several times the initial one and represent a large 
fraction of the final stiffness of the cemented sample. 
When samples are previously compacted, the initial host- 
ing medium for the cement becomes stiffer (despite frac- 
turing). In this case, further addition of small amounts 
of cement can represent an smaller fraction of the final 
stiffness of the sample. 

Results in Figure |H1 and Figure El suggest that contact 
cement theory may work well for samples which are ce- 
mented after some compaction is achieved. However, it 
may overestimate sound velocity for the case of cement 
added to weakly compacted media. The latter case, can 
be that of artificially cemented granular materials such 
as young sediments cemented with perforation mud, or 
laboratory produced materials. We believe that the case 
where cement is added after some compaction is the most 
realistic for rocks, since chemical cementation is achieved 
after some pressure-temperature conditions are reached 
with burial. 

This behaviour is further illustrated in FigurelTUlwhere 
results are shown for several cases in which cementation 
starts after some compaction. In this Figure, each trend 
may represent a possible diagenetic path followed by a 
granular material. The starting framework in all cases is 
that of a very loose material of porosity <fi = 41.2%, rep- 
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FIG. 9: Pulse velocity for different cement saturations. Un- 
cemented sample is compacted with SOMPa confining pres- 
sure. After compaction, cement is added preferentially near 
contacts. Theoretical predictions can be good estimates of 
velocity bounds for precompacted samples. 
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FIG. 10: Each curve represents a possible diagenetic path 
for a granular material. Samples are first compacted to a 
porosity (j> due to burial. Then, cementation process begin. 
The more compacted the samples, the smaller the relative 
increase of velocity with S c - 



resenting sediments after settling. The initial loose pack 
is then compacted to a porosity r/> , simulating the me- 
chanical effects of burial. After compaction, a relatively 
hard cementing material is added in varying amounts. 
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FIG. 11: Velocity increase in soft samples depends mainly 
in the amount of cement. For stiff samples, velocity increase 
depends strongly on the elastic properties of cement. 

As shown in Figure 1111 we observed the same quali- 
tative behaviour when considering cementing materials 
with different elastic properties. In the Figure, sound ve- 
locity in the uncompacted sample increases by a factor 
of w 10 when a relatively hard cement saturation varies 
from — 0.1. For the same cement saturation, a relatively 
soft cement increases sound velocity in a factor of i=s 7. 
In the case of a compacted (stiffer) sample, this variation 
approximates to ~ 25% in the case of the hard cement 
and « 12% for the soft cement. According with these re- 
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suits, velocity rise in soft samples mainly depends on the 
amount of cement while in stiff samples the most relevant 
factor is the type of cement. 
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FIG. 12: Acoustical response depends on microstructure 
details. Contact cement leads to a downward concavity in 
Vp(Sc) trend with decreasing slope as cement saturation in- 
creases. Friable cement leads to the opposite behaviour. 
Coating cement corresponds to an intermediate case. 

Cements can also accumulate as a coating cement on 
grain surfaces, or at pore bodies far from contacts. The 
former, can be the case when cement precipitates as a 
solid material in low flux zones (Manmath and Lake, 
1995). The latter case, is encountered in some sands 
which held together by confining pressure only. Following 
(Dvorkin and Nur, 1996), this kind of cement is termed 
here as friable or pore body cement. This geometric char- 
acter of cementation process can be explored in our simu- 
lations by following details given in section lTlIDI for coat- 
ing cement and section IIII El for friable cement. Methods 
to simulate the case of cement at contacts, treated pre- 
viously, are discused in section ITlI CI 

In Figure^] the lowest velocity trend is obtained when 
solid cement is added preferentially as pore-filling ma- 
terial. Since this kind of cement accumulates far from 
contacts, it contributes marginally to rock stiffness. The 
main effect of this kind of cement is to reduce porosity. 
As more pore-filling material is added, cement particles 
start to bridge separated grains. Then, blocks of cement 
particles are also compressed during pulse propagation 
with the net effect of an increasing slope of velocity trend 
as porosity reduces. 

It is to be noticed that if clays were added during 
settling, we should expect a different behaviour. In 
that case, above a critical concentration of clays (about 
« 40%), they become load-bearing and completely sur- 
round disconnected quartz grains. In such case, soft clays 
dominate the acoustics and velocity decrease for porosi- 
ties below some thereshold value (Hezhu, 1992). 

In the case treated here, Figure IT2"1 shows a monotonic 
velocity increase with cement saturation. Here, clays 



or any other cementing material is added after settling. 
Therefore, quartz grains remain in contact in the whole 
porosity range. As a consequence, soft clays may reduce 
porosity but sound velocity can only decrease slightly 
due to inertial effects. This result suggests that for low 
porosity sands with high clay content, high sound veloc- 
ity should be expected if clays are of diagenetical origin. 
Low porosity-velocity can be related to the depositional 
origin of clays. This may be an important factor to con- 
sider when reconstructing the history of a reservoir sam- 
ple. 

The velocity trend for coating cement lies in between 
limiting cases of friable and contact cement. In this 
fraction of cement deposits near contacts and 
contributes to rock stiffness, while a significant volume 
of cement mainly reduces porosity. Figure IT2l shows that 
for the same porosity, this kind of cementation leads to 
higher velocities than those of friable sand and lower 
velocities than those of the contact-cemented sample. 
These results show that microstructural details are very 
important parameters which influence in the velocity- 
porosity trends of cemented samples. 

To further compare our results with available experi- 
mental data, we calculated the P-Wave modulus from re- 
sults in Figure Computed moduli are shown against 
p " ' □ 
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FIG. 13: In the medium to high cement saturation, contact 
cement leads to a quasilinear trend for P-Wave modulus with 
porosity in agreement with experimental results (Dvorkin et 
al., 1994b). Coating cement and friable sands shows an up- 
ward concavity observed in some naturally occurring sands 
(Dvorkin and Nur, 1996). 

As expected, results in FigurefHflshow different trends 
for P-Wave modulus depending on the amount of ce- 
ment and cementation scheme. In the case of contact 
cement, a quasilinear trend is observed for cement satu- 
rations above » 10%. This result, qualitatively coincides 
with previous experimental results reported in (Dvorkin 
et al., 1994b) and references therein. Similar behavior 
was observed for the case of coating cement. Neverthe- 
less, P-modulus is smaller for the range of cement sat- 



uration studied. Velocity trends obtained for the case 
of friable cement can be compared with experimental 
data obtained for some samples with scarce intergran- 
ular cement (Dvorkin and Nur, 1996). These results are 
in good agreement with the varios correlations between 
P-modulus and porosity observed in reservoir sands. 

It must be noted that details of microstructure are 
hardly included in most theoretical approximations. 
Some models treat the case of two a phase composite 
and/or spherical inclusions, being useless in the range 
0% < S c < 100%. For instance, Hashin-Shtrikman model 
establish bounds to the maximum and minimum elas- 
tic moduli at zero porosity (S c = 100%). Although a 
generalization of this model would allow us to calculate 
bounds for the moduli at intermediate porosities (Mavko 
et al., 1988), in our simulations and in the experiments 
considered, voids arc filled with air with negligible elastic 
modulus. In this case, lower bound for the moduli is triv- 
ial(zero) in all the range of cement saturation considered. 
Although Contact Cement Model explicitly includes the 
microstructure, some the assumptions made in the the- 
ory are not supposed to be valid for large amounts of 
cement. 



V. CONCLUSIONS 

Simulation results show that cements modify elastic 
properties of sands in different manner, depending on 
the relative stiffness of cement and the hosting medium, 
the amount of cement and its localisation within pores. 

If cement is added preferentially near contacts, it can 
contribute noticeably to the stiffness of the composite, 
the effect being enhanced in uncompacted granular sam- 
ples. Friable cement leads to relatively low velocity sam- 
ples even when for moderate cement saturations. Veloc- 
ity trends show characteristic concavities in each case. 

As confining pressure increases, the hosting medium 
for the cement becomes stiffer by compaction and ce- 
ments have an smaller effect on the sound velocity 
of the final composite. As a result, velocity-porosity 
trends differ for samples that follow different compaction- 
cementation paths. Although contact cement theory does 
not account for the confining pressure effects, it can give 
a good estimate for the sound velocity in precompacted 
samples with small amounts of contact cement(5' c < 
15%). 

Simulation results presented, qualitatively and quanti- 
tatively reproduce the available experimental data where 
porosity reduces and velocity increases due to cemen- 
tation. The proposed methods, are suited to take into 
account the microstructural details of cementation pro- 
cess and their acoustical implications. As these details 
are hardly included in most effective medium approxi- 
mations, simulation techniques can be considered an al- 
ternative tool to model the acoustics of cemented sands 
and to get new insights on the underlying physics behind 
the phenomena. 
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Nomenclature 

FZI Geometric index to characterise microstructure 
P s Pore shape related value 
S v , Pore specific surface 
t Tortuosity 
F c Contact force 

F„ Hertz Force for a pair of grains in contact 

F s Shear force of a grain-to-grain contact 

AF S Change in Shear force due to a perturbation 

AF„ Change in Normal force due to a perturbation 

AF C = AF n + AF„ 

a c Radius of the contact area 

Ei,Ef Young modulus of grain i and effective Young 
modulus for a contact 

K n Hertz normal stiffness of a contact 

n s Shear stiffness of a grain-to-grain contact 

fii2 Unitary vector joining the centre of grains 1, 2. 

s Unitary vector tangential to the contact 

v,Vi,v c Poisson coefficient of grains, particle i and ce- 
ment 

i?, Rf Grain radius and effective radius for Hertzian con- 
tacts 

r. Position vector of grain i 

(i Poisson coefficient of dry friction between grains 

e n Restitution coefficient for normal interaction 

7„,7s Damping constants for normal and tangential di- 
rections to the contact. 

c w Width of cells that discretise the sample 

Pa Pg Volumetric mass density of cement and grains 

£ Overlapping between spheres in contact 

£o Overlapping between spheres in mechanical equilib- 
rium 

£ Relative displacement of contacting particles in direc- 
tion tangential to the contact 

Co Relative displacement of contacting particles in direc- 
tion tangential to the contact when particles are in 
mechanical equilibrium 
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G, G c Shear modulus of material composing grains and 
cement 

A c Lame constant of cementing material 

C n Normal stiffness for cement particles interaction 

C Elastic constant for the interaction between cement 
particles 

X12 Vector joining particles in the undisturbed cement 
lattice 

U12 Relative displacement of particles from their posi- 
tions on the undisturbed lattice. 

m c Mass of a cement particle 

Vp C , V sc P-wave and S-Wave velocities of cementing ma- 
terial 

V p P-wave velocity of the grains and cement composite 

P c Confining pressure 

4> Porosity of the cemented sample. 

4> Porosity before cementation. 

S c . Cement saturation 

Kf,Gf, M c Bulk, Shear and Compressional modulus of 
cement 

£1 Computational cost 

t c Characteristic time for the interaction among parti- 
cles 
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APPENDIX A: EFFECTIVE MEDIUM 
APPROXIMATION FOR CEMENTED SANDS 

According to the theoretical model proposed by 
Dvorkin et al. (1994a); Dvorkin et al. (1994b) and 
Dvorkin and Nur(1996), the initial porosity <j) of an un- 
cemented sample is decreased to cj> by the addition of 
cementing material. Once a given volume of cement is 
added, the effective bulk Kf and shear Gf moduli of the 
composite is calculated according with Eq. IA1I 



where p c ,V pc ,V sc are mass density, P-wave and S-wave 
velocities of cementing material and N is the average 
number of contacts per grain. The normal and tangen- 
tial stiffness of the contacts in the cemented sample, de- 
noted as S n and St respectively, are given as function of 
porosity change due to cementation in Eq. IA2I below. 



S n 
St 



■- A n a 2 
Am 2 - 



- B n a 



C s 



(A2) 



where: 



/ 2 4, -<p \ 

\3 1-00 J 



9 I l 4>o-4> 
A 1 3JV l-0o 



1/2 



1/4 



for cement in layers 



for cement at contacts 



(A3) 

The coefficients of the a powers in Eq. IA2I are given 
in formulae IA4I 



A n = -0.024153A- 1 - 3646 
B n = 0.20405A- 089008 
C n = 0.00024649A- 19864 



A s 
B s 

C s =- 
where 



-10~ 2 (2.26i/ 2 + 2.07^ + 2.3)A^ 
(0.0573z/ 2 + 0.0937i/ + 0.202)A^ 2 
-10~ 4 (9.654i/ 2 + 4.945^ + 3.1)A^ 3 



(A4) 



tJi = (0.08167^ 2 + 0.4011i/ - 1.8186) 
uj 2 = (0.027^ 2 + 0.0529j/ - 0.8765) 
w 3 = (0.079zy 2 + 0.1754z/ - 1.342) 



A 



ttG 
A - Gc 
JVs ttG' 



2G C {l-v){l-u c ) 



l-2v c 



where v, v c are respectively the grains and cement Pois- 
son ratio, G is the shear modulus of grains and the rest 
of symbols have the same denotation than in previous 
equation. 

The theoretical predictions for the P-Wave velocity are 
then obtained for the cemented sample by substituting 
elastic moduli given in Eq. lAll in linear elasticity Eq. IA5I 



if / + 4/3G / y /2 



Pf 



(A5) 



where V p is the sound velocity of the composite system of 
grains and cement and pf the effective mass density. In 
our simulations, we have used Eq. IA5l to fit the obtained 
data. 



K f - |JV(1 - <j> )M c S n 

Gf = lK } + ^N{l-4> )G c S t 



M c = V pc p c 



G c — Vsc Pci 



(Al) 



APPENDIX B: THE EFFECT OF RESOLUTION 
IN THE COMPUTATIONAL COST 

Resolution Q, affects computational cost in three dif- 
ferent ways. First, cells width c w w i?/S7, being R the 
average grain radii. As £1 increases, c w reduces and the 
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number of cells to discretise the sample increases as w ft 3 
in three dimensions. The larger the number of cells the 
larger the memory needed to store the information re- 
lated with each cell. Additionally, as the number of ce- 
ment particles raise, the number of operations to com- 
pute forces, velocities and actualise positions increases 
proportionally. 

On the other hand, for smaller cells widths c w , the 
mass of cement particles goes as m c rts f2~ 3 and elastic 



stiffness constants goes as k w O -1 . The characteris- 
tic integration step t c is determined by the characteristic 
time related to interactions t c w (m c / ' n) 1 / 2 rts 0. There- 
fore, as resolution ft increases, the number of integration 
steps to simulate pulse propagation in the system dur- 
ing a unitary time rises linearly. As both the number 
of cement particles and integration steps increase, the 
computational cost is higher. 
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